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ABSTRACT 

We compute trajectories of dust grains starting from a homogeneous surface activity- 
profile on a irregularly shaped cometary nucleus. Despite the initially homogeneous dust 
distribution a collimation in jet-like structures becomes visible. The fine structure is 
caused by concave topographical features with similar bundles of normal vectors. The 
model incorporates accurately determined gravitational forces, rotation of the nucleus, 
and gas-dust interaction. Jet-like dust structures are obtained for a wide range of gas- 
dust interactions. For the comet 67P/Churyumov-Gerasimenko, we derive the global 
dust distribution around the nucleus and find several areas of agreement between the 
homogeneous dust emission model and the Rosetta observation of dust jets, including 
velocity-dependent bending of trajectories. 


Introduction 


The modeling of the dust and gas emissions of cometary nuclei requires linking physical and 
chemical processes ranging from the formation and release of molecules to the emission of macro¬ 
scopic particles, driven by solar radiation and gravitational forces. Advanced models relating gas 
emission and subsequent acceleration of dust grains from the surface layer of a comet have been 


developed over the years (Zakharov et al. (2009); Combi et al. (2012)), but the detailed mechanism 
for dust-jet formation is not well understood |Vincent et ah] (2015). Space-craft observations start¬ 
ing with the Halley flyby 1986 reveal complex jet/filament structures near the cometary nucleus 
(Keller et al. ( 1994[ )). The Rosetta mission ( |Schnlz| ( |2009[ )) observations cover several months of 
imaging the activity of comet 67P/Churyumov-Gerasimenko (67P/C-G) and have started to release 
images captured with its navigation cameras (NAVCAM), revealing dust jets originating from the 
surface (Rosetta blog ( 2015D ). Compact dust particles originating from 67P/C-G with sizes from 
0.03 to 1 mm and fluffy dust aggregates (sizes 0.2 to 2.5 mm) have been registered with the Grain 
Impact Analyzer (GIADA) instrument on ESA/Rosetta (Fulle et al. ( 2015| )) and give typical dust 
velocities in the range of 2 to 6 m/s from August 2014 to January 2015 at the Rosetta distance 
from the nucleus 25 — 50 km). 
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The collimation of dust particles into jets has been previously linked to surface structures by 
either introducing isolated areas of cometary activity with respect to gas emission (“spots”) ( [Combi 


et al. ( 2012| )), topographical features such as pits and floors of craters (Keller et al. ( |1994| )), 


or 


terraced regions ( jFarnham et al. ( 2013| )). 


In Sect. we propose a different collimation mechanism based on the assumption of a ho¬ 
mogeneous surface activity with dust emitted along the surface normal. The model puts strong 
constraints on the origin of the dust particles. We obtain the global dust distribution within sev¬ 
eral cometary radii without introducing additional parameters or localized sources of gas and dust 
emission. A detailed topographical model of the comet is required to account for the concave ter¬ 
rain. The collimation and formation of dust jets in the homogeneous activity model result from the 
convergence of trajectories originating from concave areas. In addition, comparison with Rosetta 
images reveals the importance of cometary rotation, which causes a velocity dependent bending 
of jets due to the Coriolis effect. The Coriolis effect allows one to infer the absolute dust velocity 
and to rule out scenarios of dust emitted with zero initial velocity from the surface. The tenuous 
gas atmosphere around the comet further accelerates emitted dust, but is not the cause of the jet 
formation. The dust jets are robust over variations of four orders of magnitude of the gas-dust 
interaction. 


We start by modeling the gas density and velocity field around the comet within a simple 
model satisfying the equation of continuity for an expanding gas cloud. In the second step, we 
consider the emission of homogeneously distributed dust grains from the cometary surface and 
demonstrate a mechanism for collimated dust jets despite uniform emission. In contrast to previous 
models of cometary jets, we fully incorporate the rotation and shape of 67/C-P by working with 
a detailed topographical model of the nucleus. The surface topography is directly imprinted in 
the dust distribution around the nucleus and comparison of model with Rosetta images allows one 
to conclude that dust is leaving the surface with a non-zero velocity already acquired within the 
outer mantle. It is then directed on average along the surface normal. In addition the Coriolis force 
results in the bending of jets and puts strong constraints on dust velocities to be consistent with the 
observations of curved jets from 67P/C-G. In Sect. we incorporate insolation conditions on the 
surface and take the position of the Rosetta spacecraft into account to compare column densities 
obtained from the model with Rosetta NAVCAM images. Jet-like structures are reproduced at 
locations seen in the NAVCAM images and the comparison of model and observation allows one 
to determine the dust velocity in addition to verify the proposed collimation mechanism. The 
model itself is evaluated efficiently on massively parallel processors; a complete implementation is 
provided online Q 


^http://github.org/noma/covis 
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2. Homogeneous surface-activity model of the nucleus 


Modelling radiation-driven activity of comets requires setting up a three-dimensional model 
of the physical (temperature, density/porosity, gas/solid phase) and chemical properties (molecu¬ 
lar reactions) throughout the nucleus and coma (Huebner et al. ( |2006 )). In computational fluid 
dynamics approaches, the nucleus is subdivided into volumetric cells and the diffusion and heat 


equations are solved (Huebner et al. (2006); Bruck Syal et al. (2013)). The dust is subsequently 
embedded into the gas-model and accelerated by gas drag within the porous mantel and outside 
the cometary body. The gas eventually disperses in lateral directions and the gas density falls 
of as 1/r^ with increasing distance r, while solid dust grains remain compact particles under the 
influence of external forces. 

For the homogeneous activity model, we start by computing the gas density and gas velocity 
field around the nucleus and subsequently consider the dust motion under gas drag. 


2.1. Gas model 


Since cometary activity is largely thermally driven, differences in gas and dust production 
between the night and day sides of the comet are expected. Our aim is to model the dust distribu¬ 
tion, which is linked to the gas release inside the comet. In a first approximation, we consider the 
day/night hemispheres as completely active/inactive areas for the dust. For simplicity, we consider 
gas emanating around the whole comet. To obtain a homogeneous sampling of the comet with 
respect to the number of emitted particles (either dust or gas molecules) we employ a surface mesh 
with equal-area mesh-cells. This approach has advantages for an efficient parallelization of the 
required computations. 


Our starting point is the triangular surface mesh of the comet with detailed three-dimensional 
topography, including concave areas. Since no shape model of 67P/C-G has been published yet. 


we use a remeshed model based on the higher resolution model obtained by M. Maimer (Maimer 


( 2015| )) based on Rosetta-NAVCAM pictures. The shape model is uniformly remeshed using the 
Surface Remesh module ( [Zilske et al. (2008); Surazhsky & Gotsman ( 2003[ )) of the Amira software 
(Stalling et al. (2004)). The resulting mesh in Fig. left panel, contains 19806 triangles with a 
mean area of 2379 m^. The key property of the new mesh is the small spread of the areas of the 
mesh triangles, as shown in the histogram (Fig. right panel). This allows one to simulate a 
homogeneous surface emission rate by assigning a fixed number of representative test particles to 
each cell with well-defined surface normal direction and uniform areal coverage. In our simulation 
we set this number to 1. Surface location specific activity profiles depending on local sunlight 
conditions or active regions are optionally included by weighting the test particles accordingly, 
without requiring a new computation of trajectories. 


In a first calculation we consider the test particles to be representatives of gas parcels emitted 
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from a porous surface. The number density of the gas A^gas 


is considered later for five different 
cases ranging from lO^^-lO^^ m“^ at the surface, corresponding to low and high cometary activity 
of 67P/C-G (Combi et al. ( 2012| )). The gas velocity Vgsis 
et al. ( |2006 ), eq. 3.33) 


is taken as the most-likely speed (Huebner 


^th = 


'kBT 


m. 


( 1 ) 


gas 


resulting in v^h ~ 150 m/s for CO 2 molecules at T = 130 K, a velocity considered by Combi et al 


(2012) on a spherical, non-rotating nucleus. The dust and gas distribution does not depend on 


the precise values of these parameters, as shown later by varying the gas-dust interaction by four 
orders of magnitude. The initial direction of the velocity 'L’init of the gas emanated at distance /iinit 
from the center of the mesh cell is taken along the outward surface normal Hi of the mesh-cell i: 


Vi{t — 0) — '^init '^2 5 — 1? • • • 7 -^cells 

ri{t = 0) = (center of triangle(i)) -h /iinit 


( 2 ) 

(3) 


The positions and velocities of the gas-parcels are specified in the body-fixed (i.e. rotating) frame. 
The comet rotates with period T^ot in the inertial frame around the e^-axis (see Fig.[^. Each parcel 
is numerically integrated in the body-fixed frame according to the acceleration 


^gas 


(f) = 


(jJ = 


^grav T nc0ntrifugal T nQoriolis 

—V(/)(r) — uj X (uj X r) — 2uj X V 

(0,0, 27r/Trot), 


(4) 

(5) 

( 6 ) 


where agrav denotes the contribution due to gravity computed from the polyhedral shape model 


following Conway (2014), and acentrigaU ncorioiis account for the body-fixed frame of integration. The 
lateral expansion of gas is accounted for by a subsequent Caussian convolution, and the mixing of 
gas stemming from different regions of the concave object by a binning and averaging of the velocity 
direction of ^ 10^ position and velocities. The positions and momenta are obtained by integrating 
Eq. The binning is performed with a cubic mesh with cell length 0.2 km extending ±20 km 
around the center of the nucleus. This simplified atmospheric model reproduces on average the 
(r) ^ 1/r^ density decay with increasing distance r expected in the cometary coma (Edberg 


N, 


gasv 


mass 

1.0 X 10^3 

kg 

volume 

19.34 

km^ 

density 

517 

kg/m' 

rotation period 

12.4043 

h 

rotation axis right ascension a 

69.3 

0 

rotation axis declination 6 

64.1 

0 


Table 1: Properties of the model for 67P/C-G. 




















et al. ( 2015| )) and provides an upper limit for the influence of gas drag on the dust outside the 


nucleus. Finally, the gas number-density distributions are linearly scaled up to have maximal values 
^gas^ — 10^^, 10^^, 10^^, 10^^, 10^^ m“^ corresponding to five different scenarios. The density in the 
z = 0 plane and the velocity field is shown in Fig. |^for the gas velocity 150 m/s. 


2.2. Dust dynamics 


In the next computation step, we obtain the dust dynamics within the expanding gas cloud. 
The dust trajectories are obtained by integrating 

«dust(?^ = ^gas-drag T ^graiv T d!/c0ntrifugal T d!/Qoriolis (^) 

= ^Q(aA^gas(f0mgas(f7dust - f4as)|'^dust “ f^asl “ V(/)(r) - 2cd X Fdust - Cd X (cd X f), 

which includes the acceleration of dust embedded in the more rapidly expanding gas. The gas-drag 
acceleration depends on the dust-particle mass mdust ^md radius i?dust 5 well as on the momentum 
of the gas molecules. If we take = 2 (Keller et al. (1994)) and insert the gas velocity from Eq. 0 
we obtain an overall factor Cgas-drag foi* fhe gas-dust interaction 


d^-gas-drag 


= c, 


'gas-drag 


('ddust 


^gas 


) I 'ddust 


^gas I 


^gas I 


Cgas-drag — 


gas 


^dust 

^dust 


( 8 ) 

(9) 


In the following we vary Cgas-drag by five orders of magnitude by altering the gas density in the 
range 10^^ — 10^^ m“^ for dust particles of size 1 /im and dust density 1000 kg/m^. Equivalently, 
the parametric study corresponds to changes in the gas temperature, or/and the particle size and 
mass over an equivalently wide range. 

In addition to the gas drag encountered by dust outside the nucleus, different scenarios for 
acceleration of dust already in the outer mantle of the nucleus have been considered by [Hnebner 


et al. (2006), chapter 3.4, who model dust flow and acceleration in the pores of the mantle. In this 


case, dust is ejected out of the nucleus with finite initial velocity, gained within less than a meter. 

One direct consequence of the final initial velocity model is that the local topography of the 
cometary surface is reflected in the initial velocity direction, which is then expected to be (on 
average) along the surface normal. If on the other hand the initial dust velocity on the surface is 


taken to be zero as suggested by Combi et al. (2012), no surface normal direction is imprinted in 
the dust field. 


As we demonstrate next, only the finite initial-velocity model leads to jet formation from a 
homogeneous surface-emission of dust. The homogenous dust emission is implemented by assigning 
one dust particle to each center of the equally-sized triangles of the surface mesh. Three different 
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dust-velocities are considered, |'r’dust,init I = 0,1,2 m/s. The direction of the non-zero velocity cases 
is pointing along the surface normal of each triangle. 


3. Results for 67P/Churyumov-Gerasimenko 


It is instructive to study the relative contribution of the different terms contained in Eq. 0 
for the specific case of 67P/C-G. We consider as exemplary particle a dust grain emitted from 
the almost fiat surface of the big lobe of the comet moving in the gas field = 10^^ m“^, 

Vg^i = 150 m/s) shown in Fig. The velocity and position of the dust particle are determined 
by the gas drag, the Coriolis effect, gravitational and centrifugal forces, shown in Fig. along the 
trajectory. While the infiuence of gas drag drops rapidly with increasing distance, the rotation of the 
comet has a large infiuence on the dust trajectories through the velocity-dependent Coriolis effect. 
For 67P/C-G the rotation of the nucleus cannot be neglected as was done in previous studies in 
view of the small relative value of the centrifugal force (Zakharov et al. ( 2009| )). If only the Coriolis 
force is present, the motion of a particle with its velocity vector pointing radially outward in the 
plane orthogonal to the rotation axis follows a circle with radius 

I'^dust ) I 


-f^Coriolis — 


2u 


( 10 ) 


The NAVCAM observations allow one to determine i?Corioiis ^.nd to validate the dust model dis¬ 


cussed here, see Sect. 3.3 


3.1. Sunlight driven activity model 


A direct comparison with Rosetta observations requires considering the solar illumination 
condition and spacecraft position with respect to the comet. We construct a directional sunlight 
model using the sun-comet vector obtained from the HORIZON ephemeris file of the heliocentric 
position of the comet. The shape model is then rotated to align the rotation axis with the pole 
location in Tab. [T} Any precession of the rotation axis or change in rotation period due to torque 
acting on the comet is neglected. For the given day of the observation we rotate the shape model 
around its rotation axis to match the light conditions of the NAVCAM image frames published in 


Rosetta blog (2015), for which no exact time stamp information is given. The published NAVCAM 


images have been processed by ESA to enhance the visibility of dust and jet structures. For a direct 
comparison, we compute the comet shadow on itself and the shadow on emitted dust particles. In 
addition, we use a sunlight-driven activity model for dust emission and remove all dust trajectories 
that originate from shadowed surface cells. This simplified model neglects border effects at the 
day/night terminator and thermal inertia, leading to a prolonged activity in previously illuminated 
areas. For the demonstration of the collimation effect, this simplified activity model suffices in the 
vicinity of the nucleus where dust spends about 1 h (see Fig. [^, compared to the comet ary rotation 
period of 12 h. 
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3.2. Dust column densities 

To compare the dust model with Rosetta NAVCAM observations of 67P/C-G we obtain the 
dust column density along the camera direction extending ±20 km from the center of the nucleus. 
We use two equivalent methods to show the column density: in Fig. left column, we display 
the binned dust density in a grid with cell sizes 0.1 kmx0.1 kmx40 km, where the 40 km axis is 
oriented along the NAVCAM view direction, and show for comparison the dust density along a 
slice through the center of the comet. Alternatively, a linear superposition of point clouds taken at 
equidistant time-intervals and represented by finite-sized pixels yields the column density shown in 
Figs. [^and[^ The latter method provides a phase-space portrait of the velocity field around the 
nucleus at various positions. 


3.3. Coriolis effect 


Curved trajectories are clearly discernible in the NAVCAM images from May 3 and 12, 2015, 
see Fig. A comparison with computed dust trajectories leads to an estimate of the velocity 
3 — 4 m/s by the curvature of the trajectory (Eq. in agreement with the velocities seen along 

the trajectories for the cases of initial dust velocities in the range of 1 — 2 m/s (right column of 
Fig.g. The velocities are in agreement with CIADA measurements at r ^ 25 km reported by 


Fulle et al.| ( |2015[ ). The Coriolis effect puts strong constraints on the initial dust velocity. For 
the zero-velocity case the slow speed of the dust at lift-off leads to a hovering of dust above the 
surface and a large sidewards drift due to the Coriolis force. With an increase in gas-density a 
faster velocity gain is possible, but the dust then follows essentially the gas distribution and the 
curvature and collimation effect discussed next are both absent (Fig. lower right panel). 


4. Results 

The computed results show fine-structured and collimated jets within the dust distribution 
and the column dust density. The dust jets are forming despite the homogeneous dust release 
across the entire (sunlit) surface and are naturally arising from the simulation, without need to 
assign isolated patches of activity. In the here-demonstrated collimation mechanism a homogeneous 
emission of dust particles from the surface with finite velocities along the local surface normal is the 
pre-requisite for jet formation. The collimation is driven by the cometary surface topography and 
reflects concave (crater-like and larger scale) surface features, whose surface normals intersect in a 
focal area several km away from the surface. The collimation disappears if the gas velocity field 
acts as a dispersive, smoothening agent. For the gas-densities studied, collimation is lost at gas 
densities exceeding > 10^^ m“^ (or equivalent changes of the other parameters in Eq. Q), 

see the sequence covering five orders of gas-dust interactions in Fig. 





The collimation and bundling of trajectories, which are initially more uniformly distributed 
either in position or momentum space, is a generic physical feature seen whenever a perturbation 
slightly distorts an otherwise uniformly propagating trajectory field. The phase-space caustics due 
to concave terrain and the deflection of particle flux by the Coriolis force lead to the piling up 
of particle density into jets. The point cloud shown at equidistant time intervals in Fig. (upper 
right panel) reveals the fronts of the dust particles with converging trajectories. Similar effects have 
been observed on the nanometer scale leading to branched electron flow over a random potential 
(Topinka et al. ( |2001[ )) and on the km scale for freak waves in otherwise almost Gaussian seas 
dHeller et al] ( [2008| )). 


5. Conclusions 


We have shown that homogeneous gas and dust emission on a detailed topographic model of 
the rotating cometary nucleus accounts for many jet-like features seen in NAVCAM observations 
of 67P/C-G. In contrast to other models discussed, we accurately compute the gravity around 
concave regions (including the neck region of 67P/C-G) and consider the rotation of the nucleus. 
We relate jet-like collimated dust densities up to several km away from the cometary surface to 
concave regions of the general terrain, but without assuming a higher per-surface activity within an 


isolated region inside or around the crater as done by Keller et al. (1994) or localized gas activity 


on an otherwise smooth cometary surface (Combi et al. ( 2012| )). The collimation into jet-like 
dust structures is purely caused by convergence of trajectories of dust emitted from cracks and 
pores in the mantle along surface normals with an initial velocity of about 1 — 2 m/s. The jets 
therefore directly reflect the surface topography. This collimation mechanism persists across four 
orders of gas-dust interaction, but disappears when dust with zero initial velocity on the cometary 
surface is considered. In addition, the Coriolis effect dominates the curvature of the trajectories 
seen in the NAVCAM images of 67P/C-G and allows one to determine the dust velocity within 
specific jets (3 — 4 m/s). Isolated regions of reduced activity and day/night effects can be readily 
incorporated in the model by removing dust trajectories originating from the corresponding surface 
cells. The algorithm for trajectory computation and visualization runs efficiently on massively 
parallel processors and is available as open access code. It could be incorporated into other gas 


and dust models, that focus on additional spatially/temporally localized dust-jet formation Vincent 


et al. (2015). 


Acknowledgements 

We thank M. Maimer for helpful discussions and providing the shape model of comet 67P/C-G 
and the ESA/Rosetta NAVCAM team for releasing images of 67P/C-G under the CC BY-SA ICO 
3.0 license on a continuous basis. TK acknowledges support by a Heisenberg fellowship of the DFG 
grant KR 2889/5 and helpful discussion with P. Kramer. Computational resources and support by 
the North-German Supercomputing Alliance (HLRN) are gratefully acknowledged. 





























- 9 - 


REFERENCES 

Bruck Syal, M., Schultz, P. H., Sunshine, J. M., et al. 2013, Icarus, 222, 610 

Combi, M. R., Tenishev, V. M., Rubin, M., Fougere, N., h Gombosi, T. 1. 2012, The Astrophysical 
Journal, 758, 144 

Conway, J. T. 2014, Celestial Mechanics and Dynamical Astronomy, 121, 17 

Edberg, N. J. T., Eriksson, A. L, Odelstad, E., et al. 2015, Geophysical Research Letters, n/a 

Farnham, T. L., Bodewits, D., Li, J. Y., et al. 2013, Icarus, 222, 540 

Fulle, M., Corte, V. D., Rotundi, A., et al. 2015, The Astrophysical Journal, 802, L12 

Heller, E. J., Kaplan, L., & Dahlen, A. 2008, Journal of Geophysical Research: Oceans, 113, 1 

Huebner, W., Benkhoff, J., Capria, M.-T., et al. 2006, Heat and Gas Diffusion in Comet Nuclei 
(Bern, Switzerland: The International Space Science Institute) 

Keller, H. U., Knollenberg, J., & Markiewicz, W. 1994, Planetary and Space Science, 42, 367 

Maimer, M. 2015, http://mattias.malmer.nu/2015/02/shapemodel-update-2/ (retrieved 5/2015) 

Rosetta blog, . . 2015, http://blogs.esa.int/rosetta/ 

Schulz, R. 2009, Solar System Research, 43, 343 

Stalling, D., Westerhoff, M., & Hege, H.-C. 2004, in The Visualization Handbook, ed. C. D. Hansen 
& C. R. Johnson (Oxford, UK: Elsevier), 749 

Surazhsky, V., & Gotsman, C. 2003, in Proceedings of the 2003 Eurographics/ACM SIGGRAPH 
Symposium on Geometry Processing, Eurographics Association, 20-30 

Topinka, M. A., LeRoy, B. J., Westervelt, R. M., et al. 2001, Nature, 410, 183 

Vincent, J.-B., Bodewits, D., Besse, S., et al. 2015, Nature, 523, 63 

Zakharov, V. V., Rodionov, A. V., Lukianov, G. A., & Crifo, J. F. 2009, Icarus, 201, 358 

Zilske, M., Lamecker, H., h Zachow, S. 2008, in Eurographics 2008 Annex to the Conf. Proc., 
207-211 


This preprint was prepared with the A AS lATLX macros v5.2. 



triangle area histogram 


- 10 - 



Fig. 1.— Left panel: Triangular mesh used for the computation; color denotes the local effective 
potential (blue: 0.33 N/m, red: 0.5 N/m); the blue arrow along indicates the axis of rotation. 
Right panel: Histogram of the triangle area in the remeshed surface model (average cell area 
2379 m2). 
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Fig. 2.— Gas number density A^gas around 67P/C-G (slice at z = 0), gas velocities |Fgas| = 150 m/s, 
= 10^^ m“^. The arrows indicate the velocity field. 
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Fig. 3.— Left panel: accelerating forces acting on a dust grain emanated from the large lobe. 
Shown are the contributions due to gas drag, Coriolis effect (acting orthogonal to the velocity), 
gravitation, and centrifugal terms. Right panel: velocity and distance of the dust particle within 
the first 2500 s along the same trajectory. 


velocity (m/s) 







- 12 - 


ESA NAVCAM May 12, 2015 
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Fig. 4.— Upper left panel: ESA NAVCAM image of 67P/C-G on May 12, 2015, inverted greyscale 
to highlight curved jets (image credits ESA/Rosetta/NAVCAM CC BY-SA ICO 3.0). The cur¬ 
vature of the jets follows a circle of approximately 11.5 km (dotted arc), which corresponds to 
'^dust = 3.3 m/s. Left panels: column density from the homogeneous activity model and density 
at a plane intersecting the center of the nucleus. Right panels: Sequence of decreasing cometary 
surface dust velocities (2,1,0 m/s). The zero velocity case does not show collimated jets. 
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Fig. 5.— Upper left panel: ESA image on 20 April 2015, credits ESA/Rosetta/NAVCAM CC 
BY-SA IGO 3.0. Other panels: comparison of five different surface gas densities. Best agreement 
with theory is observed around = 10^^ — 10^^ m“^. 









